{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling 2: Create a User Defined Model using astropy.modeling\n",
    "\n",
    "## Authors\n",
    "Rocio Kiman, Lia Corrales, Zé Vinícius, Stephanie T. Douglas\n",
    "\n",
    "## Learning Goals\n",
    "* Define a new model with `astropy` \n",
    "* Identify cases were a user-defined model could be useful\n",
    "* Define models in two different ways:\n",
    "    * Compound models\n",
    "    * Custom models\n",
    "\n",
    "This tutorial assumes the student knows how to fit data using `astropy.modeling`. This topic is covered in the [Models-Quick-Fit tutorial](http://www.astropy.org/astropy-tutorials/rst-tutorials/Models-Quick-Fit.html).\n",
    "\n",
    "## Keywords\n",
    "modeling, FITS, astrostatistics, matplotlib, model fitting, error bars, scatter plots\n",
    "\n",
    "## Summary\n",
    "In this tutorial, we will learn how to define a new model in two ways: with a compound model and with a custom model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**Note:** This tutorial assumes you have already gone through \n",
    "    [Modeling 1](http://www.astropy.org/astropy-tutorials/rst-tutorials/Models-Quick-Fit.html),\n",
    "    which provides an introduction to `astropy.modeling`\n",
    "    \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "from astropy.modeling import models, fitting\n",
    "from astropy.modeling.models import custom_model\n",
    "from astropy.modeling import Fittable1DModel, Parameter\n",
    "from astroquery.sdss import SDSS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit an emission line in a stellar spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "M dwarfs are low mass stars (less than half of the mass of the sun). Currently we do not understand completely the physics inside low mass stars because they do not behave the same way higher mass stars do. For example, they stay magnetically active longer than higher mass stars. One way to measure magnetic activity is the height of the [$H\\alpha$](https://en.wikipedia.org/wiki/H-alpha) emission line. It is located at $6563$ Angstroms at the spectrum. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's search for a spectrum of an M dwarf in the Sloan Digital Sky Survey (SDSS). First, we are going to look for the spectrum in the [SDSS database](https://dr12.sdss.org/basicSpectra). SDSS has a particular way to identify the stars it observes: it uses three numbers: Plate, Fiber and MJD (Modified Julian Date). The star we are going to use has:\n",
    "* Plate: 1349\n",
    "* Fiber: 216\n",
    "* MJD: 52797\n",
    "\n",
    "So go ahead, put this numbers in the website and click on Plot to visualize the spectrum. Try to localize the $H\\alpha$ line. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could download the spectrum by hand from this website, but we are going to import it using the [SDSSClass](http://astroquery.readthedocs.io/en/latest/api/astroquery.sdss.SDSSClass.html) from [`astroquery.sdss`](https://astroquery.readthedocs.io/en/latest/sdss/sdss.html#module-astroquery.sdss). We can get the spectrum using the plate, fiber and mjd in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum = SDSS.get_spectra(plate=1349, fiberID=216, mjd=52797)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now that we have the spectrum...\n",
    "One way to check what is inside the fits file `spectrum` is the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "spectrum[1].columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To plot the spectrum we need the flux as a function of wavelength (usually called lambda or $\\lambda$). Note that the wavelength is in log scale: loglam, so we calculate $10^\\lambda$ to remove this scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux = spectrum[1].data['flux']\n",
    "lam = 10**(spectrum[1].data['loglam'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find the units for flux and wavelength, we look in `fitsfile[0].header`.\n",
    "\n",
    "FITS standard requires that the header keyword 'bunit' or 'BUNIT' contains the physical units of the array values. That's where we'll find the flux units. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Units of the flux\n",
    "units_flux = spectrum[0].header['bunit']\n",
    "print(units_flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Different sources will definite wavelength information differently, so we need to check the documentation. For example, this [SDSS tutorial](https://www.sdss.org/dr12/tutorials/quicklook/#python) tells us what header keyword to look at."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Units of the wavelegth\n",
    "units_wavelength_full = spectrum[0].header['WAT1_001']\n",
    "print(units_wavelength_full)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to select only the characters of the unit we care about: Angstroms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "units_wavelength = units_wavelength_full[36:]\n",
    "print(units_wavelength)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to plot the spectrum with all the information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(lam, flux, color='k')\n",
    "plt.xlim(6300,6700)\n",
    "plt.axvline(x=6563, linestyle='--')\n",
    "plt.xlabel('Wavelength ({})'.format(units_wavelength))\n",
    "plt.ylabel('Flux ({})'.format(units_flux))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We just plotted our spectrum! Check different ranges of wavelength to see how the full spectrum looks like in comparison to the one we saw before."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit an Emission Line with a Gaussian Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The blue dashed line marks the $H\\alpha$ emission line. We can tell this is an active star because it has a strong emission line.\n",
    "\n",
    "Now, we would like to measure the height of this line. Let's use `astropy.modeling` to fit a gaussian to the $H\\alpha$ line. We are going to initialize a gaussian model at the position of the $H\\alpha$ line. The idea is that the gaussian amplitude will tell us the height of the line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gausian_model = models.Gaussian1D(1, 6563, 10)\n",
    "fitter = fitting.LevMarLSQFitter()\n",
    "gaussian_fit = fitter(gausian_model, lam, flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,5))\n",
    "plt.plot(lam, flux, color='k')\n",
    "plt.plot(lam, gaussian_fit(lam), color='darkorange')\n",
    "plt.xlim(6300,6700)\n",
    "plt.xlabel('Wavelength (Angstroms)')\n",
    "plt.ylabel('Flux ({})'.format(units_flux))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the fit is not doing a good job. Let's print the parameters of this fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(gaussian_fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise \n",
    "\n",
    "Go back to the previous plot and try to make the fit work. Note: **Do not spend more than 10 minutes** in this exercise. A couple of ideas to try: \n",
    "* Is it not working because of the model we chose to fit? You can find more models to use [here](http://docs.astropy.org/en/stable/modeling/#module-astropy.modeling.functional_models).\n",
    "* Is it not working because of the fitter we chose?\n",
    "* Is it not working because of the range of data we are fitting? \n",
    "* Is it not working because how we are plotting the data? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compound models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One model is not enough to make this fit work. We need to combine a couple of models to make a [compound model](http://docs.astropy.org/en/stable/modeling/#compound-models) in `astropy`. The idea is that we can add, divide or multiply models that already exist in [astropy.modeling](http://docs.astropy.org/en/stable/modeling/#models-and-fitting-astropy-modeling) and fit the compound model to our data.\n",
    "\n",
    "For our problem we are going to combine the gaussian with a polynomial of degree 1 to account for the background spectrum close to the $H\\alpha$ line. Take a look at the plot we made before to convince yourself that this is the case.\n",
    "\n",
    "Now let's make our compound model!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compound_model = models.Gaussian1D(1, 6563, 10) + models.Polynomial1D(degree=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After this point, we fit the data in exactly the same way as before, except we use a compound model instead of the gaussian model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitter = fitting.LevMarLSQFitter()\n",
    "compound_fit = fitter(compound_model, lam, flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,5))\n",
    "plt.plot(lam, flux, color='k')\n",
    "plt.plot(lam, compound_fit(lam), color='darkorange')\n",
    "plt.xlim(6300,6700)\n",
    "plt.xlabel('Wavelength (Angstroms)')\n",
    "plt.ylabel('Flux ({})'.format(units_flux))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It works! Let's take a look to the fit we just made. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(compound_fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print all the parameters in a fancy way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for x,y in zip(compound_fit.param_names, compound_fit.parameters):\n",
    "    print(x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the result includes all the fit parameters from the gaussian (mean, std and amplitude) and the two coefficients from the polynomial of degree 1. So now if we want to see just the amplitude:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compound_fit.amplitude_0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Conclusions:** What was the difference between the first simple Gaussian and the compound model? The linear model that we added up to the gaussian model allowed the base of the Gaussian fit to have a slope and a background level. Normal Gaussians go to zero at $\\pm \\inf$; this one doesn't. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fixed or bounded model parameters\n",
    "The mean value of the gaussian from our previous model indicates where the $H\\alpha$ line is. In our fit result, we can tell that it is a little off from $6563$ Angstroms. One way to fix this is to fix some of the parameters of the model. In `astropy.modeling` these are called **[fixed parameters](http://docs.astropy.org/en/stable/api/astropy.modeling.Parameter.html#astropy.modeling.Parameter.fixed)**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compound_model_fixed = models.Gaussian1D(1, 6563, 10) + models.Polynomial1D(degree=1)\n",
    "compound_model_fixed.mean_0.fixed = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's use this new model with a fixed parameter to fit the data the same way we did before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitter = fitting.LevMarLSQFitter()\n",
    "compound_fit_fixed = fitter(compound_model_fixed, lam, flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,5))\n",
    "plt.plot(lam, flux, color='k')\n",
    "plt.plot(lam, compound_fit_fixed(lam), color='darkorange')\n",
    "plt.xlim(6300,6700)\n",
    "plt.xlabel('Wavelength (Angstroms)')\n",
    "plt.ylabel('Flux ({})'.format(units_flux))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(compound_fit_fixed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see in the plot that the height of the fit does not match the $H\\alpha$ line height. What happend here is that we were too strict with the mean value, so we did not get a good fit. But the mean value is where we want it! Let's loosen this condition a little. Another thing we can do is to define a [**minimum and maximum value**](http://docs.astropy.org/en/stable/api/astropy.modeling.Parameter.html#astropy.modeling.Parameter.max) for the mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compound_model_bounded = models.Gaussian1D(1, 6563, 10) + models.Polynomial1D(degree=1)\n",
    "delta = 0.5\n",
    "compound_model_bounded.mean_0.max = 6563 + delta\n",
    "compound_model_bounded.mean_0.min = 6563 - delta\n",
    "\n",
    "fitter = fitting.LevMarLSQFitter()\n",
    "compound_fit_bounded = fitter(compound_model_bounded, lam, flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,5))\n",
    "plt.plot(lam, flux, color='k')\n",
    "plt.plot(lam, compound_fit_bounded(lam), color='darkorange')\n",
    "plt.xlim(6300,6700)\n",
    "plt.xlabel('Wavelength (Angstroms)')\n",
    "plt.ylabel('Flux ({})'.format(units_flux))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(compound_fit_bounded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Better! By loosening the condition we added to the mean value, we got a better fit and the mean of the gaussian is closer to where we want it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "Modify the value of delta to change the minimum and maximum values for the mean of the gaussian. Look for:\n",
    "* The better delta so the mean is closer to the real value of the $H\\alpha$ line.\n",
    "* What is the minimum delta for which the fit is still good according to the plot?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Custom model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What should you do if you need a model that `astropy.modeling` doesn't provide? To solve that problem, Astropy has another tool called [custom model](http://docs.astropy.org/en/stable/modeling/new.html). Using this tool, we can create any model we want. \n",
    "\n",
    "We will describe two ways to create a custom model: \n",
    "* [basic](http://docs.astropy.org/en/stable/modeling/new.html#basic-custom-models)  \n",
    "* [full](http://docs.astropy.org/en/stable/modeling/new.html#a-step-by-step-definition-of-a-1-d-gaussian-model)\n",
    "\n",
    "We use the basic custom model when we need a simple function to fit and the full custom model when we need a more complex function. Let's use an example to understand each one of the custom models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic custom model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An **Exponential Model** is not provided by Astropy models. Let's see one example of basic custom model for this case. First, let's simulate a dataset that follows an exponential:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = np.linspace(0,10,100)\n",
    "\n",
    "a = 3\n",
    "b = -2\n",
    "c = 0\n",
    "y1 = a*np.exp(b*x1+c)\n",
    "y1 += np.random.normal(0., 0.2, x1.shape)\n",
    "y1_err = np.ones(x1.shape)*0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x1 , y1, yerr=y1_err, fmt='.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can define a simple custom model by specifying which parameters we want to fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@custom_model\n",
    "def exponential(x, a=1., b=1., c=1.):\n",
    "    '''\n",
    "    f(x)=a*exp(b*x + c)\n",
    "    '''\n",
    "    return a*np.exp(b*x+c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have one more available model to use in the same way we fit data with `astropy.modeling`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_model = exponential(1.,-1.,1.)  \n",
    "fitter = fitting.LevMarLSQFitter()\n",
    "exp_fit = fitter(exp_model, x1, y1, weights = 1.0/y1_err**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x1 , y1, yerr=y1_err, fmt='.')\n",
    "plt.plot(x1, exp_fit(x1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(exp_fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fit looks good in the plot. Let's check the parameters and the Reduced Chi Square value, which will give us information about the goodness of the fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_reduced_chi_square(fit, x, y, yerr, N, n_free):\n",
    "    '''\n",
    "    fit (array) values for the fit\n",
    "    x,y,yerr (arrays) data\n",
    "    N total number of points\n",
    "    n_free number of parameters we are fitting\n",
    "    '''\n",
    "    return 1.0/(N-n_free)*sum(((fit - y)/yerr)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calc_reduced_chi_square(exp_fit(x1), x1, y1, y1_err, len(x1), 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Reduced Chi Square value is close to 1. Great! This means our fit is good, and we can corroborate it by comparing the values we got for the parameters and the ones we used to simulate the data.\n",
    "\n",
    "**Note:** Fits of non-linear parameters (like in our example) are extremely dependent on initial conditions. Pay attention to the initial conditions you select."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "Modify the initial conditions of the fit and check yourself the relation between the best fit parameters and the initial conditions for the previous example. You can check it by looking at the Reduced Chi Square value: if it gets closer to 1 the fit is better and vice versa. To compare the quality of the fits you can take note of the Reduced Chi Square value you get for each initial condition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Full custom model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if we want to use a model from `astropy.modeling`, but with a different set of parameters? One example is the [Sine Model](http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Sine1D.html#astropy.modeling.functional_models.Sine1D). It has a very particular definition of the frequency and phase. Let's define a new Sine function with a full custom model. Again, first let's create a simulated dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x2 = np.linspace(0,10,100)\n",
    "a = 3\n",
    "b = 2\n",
    "c = 4\n",
    "d = 1\n",
    "y2 = a*np.sin(b*x2+c)+d\n",
    "y2 += np.random.normal(0., 0.5, x2.shape)\n",
    "y2_err = np.ones(x2.shape)*0.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x2, y2, yerr=y2_err, fmt='.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the full custom model we can easily set the **derivative of the function**, which is used by different [fitters](http://docs.astropy.org/en/stable/modeling/#id21), for example the `LevMarLSQFitter()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SineNew(Fittable1DModel):\n",
    "    a = Parameter(default=1.)\n",
    "    b = Parameter(default=1.)\n",
    "    c = Parameter(default=1.)\n",
    "    d = Parameter(default=1.)\n",
    "        \n",
    "    @staticmethod\n",
    "    def evaluate(x, a, b, c, d):\n",
    "        return a*np.sin(b*x+c)+d\n",
    "    \n",
    "    @staticmethod\n",
    "    def fit_deriv(x, a, b, c, d):\n",
    "        d_a = np.sin(b*x+c)+d\n",
    "        d_b = a*np.cos(b*x+c)*x\n",
    "        d_c = a*np.sin(b*x+c)\n",
    "        d_d = np.ones(x.shape)\n",
    "        return [d_a, d_b, d_c, d_d]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note** Defining default values for the fit parameters allows to define a model as `model=SineNew()`\n",
    "\n",
    "We are going to fit the data with our **new model**. Once more, the fit is very **sensitive to the initial conditions** due to the non-linearity of the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sine_model = SineNew(a=4.,b=2.,c=4.,d=0.)  \n",
    "fitter = fitting.LevMarLSQFitter()\n",
    "sine_fit = fitter(sine_model, x2, y2, weights = 1.0/y2_err**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x2, y2, yerr=y2_err, fmt='.')\n",
    "plt.plot(x2,sine_fit(x2))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sine_fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calc_reduced_chi_square(sine_fit(x2), x2, y2, y2_err, len(x2), 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Reduced Chi Squared value is showing the same as the plot: this fit could be improved. The Reduced Chi Squared is not close to 1 and the fit is off by small phase."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "Play with the initial values for the last fit and improve the Reduced Chi Squared value. \n",
    "\n",
    "**Note:** A fancy way of doing this would be to code a function which iterates over different initial conditions, optimizing the Reduced Chi Squared value. No need to do it here, but feel free to try."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Custom models are also useful when we want to fit an **unusual function** to our data. As an example, create a full custom model to fit the following data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x3 = np.linspace(-2,3,100)\n",
    "y3 = x3**2* np.exp(-0.5 * (x3)**3 / 2**2)\n",
    "y3 += np.random.normal(0., 0.5, x3.shape)\n",
    "y3_err = np.ones(x3.shape)*0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x3,y3,yerr=y3_err,fmt='.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
